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Abstract 

The equilibrium properties of classical self-gravitating systems in the grand canonical 
ensemble are studied by using the correspondence with an euclidean field theory with 
infrared and ultraviolet cutoffs. It is shown that the system develops a first order phase 
transition between a low and a high density regime. In addition, due to the long range 
of the gravitational potential, the system is close to criticality within each phase, with 
the exponents of mean field theory. The coexistence of a sharp first order transition and 
critical behavior can explain both the presence of voids in large regions of the universe as 
well as the self-similar density correlations in terms of self-gravity alone. 
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1 Introduction 



One of the outstanding problems of modern cosmology is the understanding of structure 
formation in the universe jT|. Visible matter at astronomical scales appears organized in 
a hierarchy of galaxies, cluster and supercluster of galaxies, which tend to be found in 
filamentary aggregates or two-dimensional sheets that encompass large regions with much 
lower density of matter and structures. These regions, called voids, occupy a large fraction 
of the observed universe, and are also organized in a hierarchy [21 [3]. The theoretical 
investigation of void formation is receiving increasingly attention [I] . 

There is general agreement that the galaxy two point correlation function is scale 
invariant (self-similar), at least at not too large scales, decaying as a power of the distance, 
1/r 7 . The exponent 7 as well as the scale of homogeneity is still controversial. Pietronero 
and coworkers found 7 ~ 1 and claim that the power law is obeyed up to the deepest 
distances $ t . Other authors gave a different exponent (7 « 1.8) and support the existence 
of an observed scale of homogeneity 00 IB]- Self-similar behaviour has also been found 
in the interstellar medium (§1 EH EEU EE1! • 

Although certainly the dynamics must be very important in order to explain these 
facts, and many physical effects might play a prominent role, it is well possible that many 
aspects of the observed structures may be understood in terms of the equilibrium states 
of self-gravitating matter alone, as claimed by the authors of Refs. |13l I14j . Indeed, it 
has been tried to apply thermodynamics to astrophysical systems since a long time (see 
for instance Ref. in this paper, we will analyze the general phase structure of a 

non-relativistic self-gravitating system at thermal equilibrium. The results are general 
and can be applied to any such system that can be considered to be in these conditions. 
We will not be concerned here with the very interesting questions of whether thermal 
equilibrium can be reached in such systems, or the way it is attained. We just assume 
that the self-gravitating system is at thermal equilibrium. 

Let us consider a system of N classical particles of mass m confined in a region of 
volume V and interacting each other via the Newtonian gravitational potential. Its hamil- 
tonian reads 



The thermodynamics of such a system is ill defined: the entropy does not exist due to the 
singularity of the potential at short distances |19| 120] . Furthermore, even if one modifies 
the potential to remove the short distance singularity, the usual thermodynamical limit 
does not exist, since the thermodynamic potentials are not extensive, due to the long 
range gravitational force. In such cases, the micro canonical specific heat can become 
negative and different statistical ensembles are not equivalent (see Refs. (2D [221 1231 Ell 
[231 1201 123 1231 E31 [30 J. In particular, the grand canonical ensemble is dominated by 
completely collapsed configurations, whatever the chemical potential. To use the ordinary 
thermodynamical tools, the potential must also be modified at long distances 1 . 

The thermodynamics of self-gravitating system has been studied since a long time 

x It has been proposed that a kind of thermodynamical limit for systems with potentials decaying as 
1/r could be taken by considering the so called dilute regime |31|. but it can be shown that this statement 
cannot hold 32 . 
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by confining a finite system on a finite box and using mean field theory. The approach 
developed in this paper is different: we will describe the thermodynamical properties of a 
self-gravitating system as the limiting case of a family of well behaved, short range systems, 
the interactions of which decay with distance more and more slowly, thus resembling more 
and more the newtonian 1/r potential. For these systems the usual thermodynamical 
limit, which is considered in this paper, does exist. To carry out this investigation, we 
will take advantage of the fact that the statistics of a self-gravitating ensemble of particles 
can be related to an euclidean field theory of a single scalar field |13M14| . in a similar way 
to the relation between the Coulomb gas and the Sine-Gordon field theory El Ev- 
idence, the remaining of the paper will rely on the techniques of euclidean field theory. 



2 Field theoretical description of a self-gravitating system 



The statistical mechanics of a self-gravitating system can be studied in the usual way if 
the gravitational potential is regularized at short and long distances. Let us choose as a 
regularized gravitational potential, denoted by U R (r), an attractive Yukawa potential, with 
range l/m-o, endowed with a hard core of size r$ at short distances: U K (r) = —Gm 2 e~ m ° r jr 
for r > ro and U K (r) = oo for r < tq. After integrating out the momenta, the grand 
canonical partition function for chemical potential \x and temperature T can be written 
as 



-'GC 



V- ( ma 2 \ 



3N/2 



\2-Kf3h 2 

where f3 = 1/ksT, a is a constant with units of length, to be specified below, and 



(2) 



-3N 



N\ 



N 



(3) 



8=1 



The partition function Zjy can be approximated by dividing the space volume in 
A = V/a 3 cells of size a (the unit length introduced above), and replacing the integrals 
by appropriate Riemann sums, which can be reorganized as a sum over cell occupation 
numbers S{, with the index i running from 1 to the number of cells, A. If a is of the order 
of the hard core size, each cell can be either void or occupied by one particle: Si = 0, 1. 
This is similar to consider self-gravitating fermions j^H]- We get 



2 n = - ex p \ -j 

5i=0 5 A =0 



U L (n-n')S n S n ,-U L (0)N 



5(£S n -N), (4) 



where b 2 = AnPGrn 2 , U L (n — re') is a proper lattice version of —U R (r — r')/(47rGm 2 ), the 
factor 1/Nl has been canceled due to the distinguishability of classical particles, and the 
Kronecker delta takes into account that the number of particles is fixed |33j . 

Since the Yukawa potential is the Green function of —V 2 + m 2 , we can choose its 
lattice version as the lattice Green function of 
the laplacian: 

U L (n - ri) 



-V 2 + ml, 
V 2 + ibq, where V 2 is a discretization of 



1 



-V 2 + rre 2 



(5) 



n,n' 
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Using © and the Hubbard-Stratonovich formula, we can write the following identity: 



b 2 

exp I — Y u ^i n ~ n ')S n S„ 



2 

n.n' 



C I [dip] exp I -i Yl ° 3 (-V 2 + rao) m - + f ' ( 6 ) 

[ n,n' n J 

where C is a number independent of S n . Plugging Eq. in (jl j) and the resulting equation 
for Zn in the expression for 2 GC , and performing the summation over Si, we can write the 
grand canonical partition function of the (regularized) self-gravitating system in terms of 
a local euclidean field theory with a single scalar field, ifjn- 

Z GC = C J[di>}ex P [-S}, (7) 

the action of which is given by 

S = \ E a 3 Vv^™„'Vv - Y ln i 1 + ^ Hn ) , (8) 

nn' n 

where g = e^e"^"^^ ^ ^"^ ) 3//2 and X> = — + tUq. Notice the unusual form of the 

action: the interaction term, ln(l + ge b ^), is unbounded from below. It behaves as —bip 
for tp — y oo. The action, however, is bounded from below if ttt-q > 0. This reflects the fact 
that the grand canonical ensemble does not exist for non-extensive (long range and/or 
purely attractive) systems. The infrared and ultraviolet cutoffs, mo and a respectively, 
make the system short ranged and repulsive at short distances (hard core of size of the 
order of a), and thermodynamics is well defined. In the naive continuum limit, in which 
a — > (g — > 0) and the relevant field configurations are smooth, we recover the action of 
Ref. 13 . This continuum action is not bounded from below and hence the corresponding 
functional integral diverges. Thus, the results of Ref. |13| are at most formal and the 
conclusions different from ours. The regulators, a and mo, play an essential role, although 
the conclusions, as we will see, are independent of them (provided they are small enough). 

The action contains four parameters: the lattice spacing, a, which obviously has dimen- 
sion of length; the inverse of the potential range, mo, which has dimension of (length) -1 ; 
b, with dimension of (length) 1 / 2 ; and g, which is dimensionless. The field is canonically 
normalized, so that its dimension is (length) -1 / 2 . It is convenient to redefine the field 
such that 4> n = bip n + ln g, and to work in terms of the dimensionless quantities k = a/b 2 , 
mg = a 2 rriQ, and h = 1/2 + /imping. Ignoring constant terms, the action reads 

<S = ^ Y ^n'Dnn'&n' ~Y ln COsh ~T ~ h ^2 ^ n ' ( 9 ) 



n iv 



where T> = — a 2 V 2 + m 2 ,. An action of the same form, with m 2 , = 0, has been obtained in 
Ref. |14j in a similar way. In that case, however, the parameter h has a different meaning. 
Notice that for h = the action © has the symmetry <p n — > —<J) n . Since the term linear 
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in (j> n breaks this symmetry, we will have ((f)) > for h > and (4>) < for h < 0. At 
h = we will have (4>) = unless the symmetry is spontaneously broken. 

The average of particle and energy densities, p and e, can be obtained as derivatives of 
the grand canonical partition function. This allows to express them in terms of the field 
(j) in the following way: 



-Y- 

V ^ 1 



+ e? 



3 1 1 
2 P ~ 2 V 



Ing-K- 1 ^ ° 



1 + e'P™ 



(10) 
(11) 



One can easily identify the first contribution to the energy in Eq. (jllj) as the kinetic 
energy and the second one as the potential energy. The later has the form ^2 n (4> n — 
Ing — K~ 1 T> nn )p n , where p n = e^ n /(l + e^ n ) is the local particle density, what implies that 
the field <f> n represents, up to an additive constant, \ng + K~ l V nn , the local gravitational 
potential 2 . 

After some algebraic manipulations, which are outlined in the appendix since they are 
not completely straightforward, we can write the following exact equation for the average 
energy, which has the expected form: 



2 P 



.-i 



■ Pn+r) 



n r^O 



(12) 



In the above equation e and p stand for the average energy and density, respectively. 
Remember that T> n n+r ~ e _m ° r /r. Notices that the fact that the term with r = is 
excluded from the sum over r in Eq. (|12|) implies that the contact term does not contribute 
to the energy. This, and the bounds < p n < 1, are manifestations of the particle hard 
core. 

3 Phase diagram 

Perturbation theory in euclidean field theory starts by identifying the minimum of the 
classical action and assuming that the relevant field configurations are small fluctuations 
around this minimum. For the action @ this is a good approximation if k is large. We 
shall argue below that it is indeed a good approximation whatever k if is small. 

The minimum of the action obviously corresponds to the constant field that minimizes 
the classical potential 

a2 



u 



In cosh — 
2 



(13) 



and, therefore, satisfies the equation 



— tanh — + h . 
2 2 



(14) 



2 Note that T> nn is independent of n due to the translational invariance of the Laplacian. 
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Figure 1: Phase diagram in the (Km,Q,g) plane. The solid line is a first order transition that ends at the 
critical point. 

The above equation has one, two or three solutions, depending on the values of Km§ 
and h. When it has one solution, it correspond to the global minimum; if there are 
two solutions, one is the global minimum and the other one a point of inflection; three 
solutions correspond to a local maximum, a local minimum, and the global minimum. In 
the last case, the minima can only be degenerate if h = 0, and then symmetry implies 
that Eq. (|14j) has either one or three solutions: for kYFTq > 1/4 the only solution is </> = 0, 
while for kiSq < 1/4 we have two solutions, = ±0 O ^ (with O positive), which are 
the two degenerate global minima, besides the symmetric solution eft = 0, which is a local 
maximum. We will denote the two global minima by (j)^ = +(f>o and (pi = —<po- The 
densities of each phase are respectively 



Hence, the solutions (f>i and 4>h describe phases of low and high density respectively. 

The phase diagram in the plane (ktSq.^), at the classical level, is displayed in figure E 
For ktuq > 1/4, the transition between the low (small g) and high (large g) density 
regimes is smooth. For ktuq < 1/4, the low density phase is separated from the high 
density phase by a first order transition which takes place at <7 PO = exp[— 1/(2ktoq)] (i.e., 
h = 0). This first order line ends at the critical point K c mg = 1/4, g c = 1/e 2 (h = 0). The 
order parameters, 8(f) = (fth — (pi and 5p = Ph — Pi, vanish at the critical point with the 
classical (mean field) exponent 1/2. This mean field critical behavior is a consequence of 
the classical approximation and could be modified by the neglected fluctuations. 

For small nm§ and h = 0, Eq. (JTIJ) gives <p h w l/(2Kfnl), a 3 p h w 1 — e~ 1 /( 2KM o) ) 
(pi ~ — l/(2Kmg), and a 3 pi ~ e -1// ( 2Km o). Thus, the low density phase is very dilute, while 
the high density phase is extremely dense. Most of the densities (including presumably 
those of physical systems) are between pi and ph, and they correspond, therefore, to the 
phase coexistence region. 



1 e*' 




(15) 



PI = 3 



a 3 1 + e^i ' 



3 1 _|_ e <t> h ' 
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4 Critical behavior and correlations 



The corrections to the classical approximation can be obtained perturbatively, by expand- 
ing the action around the corresponding minimum. For each phase we write (ft n = 4>h + '-Pn 
and (j) n = 4>i + ip n respectively, and, ignoring again constant terms, the corresponding 
actions can be written as 

S = \ E <PnV nn ><Pn> ± A ^ ip n - ^ In [l + A(e ± ^ - 1)] , (16) 

nn' n n 

where the signs + and — correspond to the low and high density phases respectively, and 
A = 1 - a 3 p h = a 3 pi. 

When ktUq is very small, A is much smaller and the actions for either phase (|16l) are 
very close to the gaussian critical point, A = 0. Therefore, the classical approximation we 
are using is very good. In addition, we have a very interesting situation: a very sharp first 
order transition that separates two phases which, in turn, are close to criticality. Hence, 
the actions (|16f) produce a very large correlation length and critical (self-similar) behavior 
over a vast range of scales. The question is to which universality class corresponds such 
critical behavior. In three dimensions, the gaussian fixed point is infrared unstable under 
perturbations of relevant operators. This means that the critical behavior of actions 
defined in the neighborhood of the gaussian fixed point can be governed by a different 
(non-gaussian) fixed point. For this to happen, the couplings of the relevant operators of 
dimension larger than one must be of the order of /otTq. In our case A <C ktuq and we are 
in the opposite case: the actions (|T6|) lie very close to the renormalized trajectory of the 
gaussian fixed point, given by A = 0, ktoq > 0. Hence, the critical behavior is governed 
by the gaussian fixed point, and thus belongs to the mean field universality class. 

The correlations can be computed to leading order in A with good accuracy. For the 
field we get (4>r4>r') — (4>r}(4>r'} = G(r — r'), with 



d 3 jp(r-r') 

a{r - r ' )= I Wfi^i- (17) 



where we have used the continuum expression since the lattice spacing is much smaller than 
the correlation length, £ = a/ {Kmfy 1 / 2 . At large distances, G(r) decays exponentially as 
exp(— |r|/£), but in a wide range of distances, a <C <C £, it is approximately self-similar, 
G(r) ~ l/\r\. 

The correlations of the density (jl()j) . T(r — r'), have the same long distance behavior 
as the correlations of the field. To leading order in A we have: 

T(r-r') = A 2 e G (°) [e G ( r - r ') - l] . (18) 

For large distances T(r) ~ G(r), and therefore the density correlations behave as l/\r\ 
for a <^ \r\ £. Hence, the present approach to self-gravitating systems predicts that 
density correlations decay as a power law with exponent 7 = 1 over a vast range of scales. 
Since £ is proportional to the assumed range of the gravitational interaction, it is well 
possible that correlations be self-similar at any observable scale. 

It is worthwhile stressing that the above analysis refers to the quasi-critical behaviour 
of each of the two phases as kWIq — ► 0. There is a true critical point at k c Mq = 1/4 and 
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g = 1/e 2 (or p c = 1/2). In this case the parameters are not small and a non-perturbative 
analysis is required in order to investigate the nature of its universality class. This critical 
point might be relevant at very high temperatures, when kWIq is of order one. 

A similar analysis can be made for small k (strong coupling regime). The symmetry 
is broken at the classical level if smg < 1/4. For small n and small m 2 , short wavelength 
fluctuations cannot induce tunneling between the two minima since it would imply huge 
fluctuations of the action, of the order 1 / kttTq . Likewise, long wavelength fluctuations 
cannot induce tunneling since the barrier is too high, of the order of 1/(8ktoq). The two 
minima describe phases of low and high density. The correlation length on each phase is 
of the order of 1/wiq. Hence, the conclusions are the same as in the weak coupling regime: 
an abrupt first order transition separates the low and high energy phases, and each phase 
present critical behaviour of mean field type. 

5 Finite volume effects: the Lane-Emden equation 

It might seem paradoxical that the mean field solution is homogeneous, since it is well 
known that the usual mean field solutions of self-gravitating systems present density pro- 
files that depend on the distance to the center of the system. There is a simple explanation: 
since we are looking for the mean field solution of a short ranged translational invariant 
system in the infinite volume limit it is natural that it be homogeneous. Had we looked 
for the minimum of the action (jHJ) on a finite box of size L we would have obtained the 
equation 

]T(-V|WVV + m|fyn - * nMn = 0- (19) 
w 

On a finite volume, the walls break translational invariance and the solution of Eq. Q19JI will 
not be homogeneous, due to the boundary conditions imposed to the operator (— V 2 L ) nn i. 
Since this operator is a discretization of the Laplacian, if L <C l/^o, Eq. (|19[) is similar 
to the isothermal Lane-Emden equation [341 135j (if mo = 0, to leading order in g we will 
have exactly the discretized Lane-Emden equation, and the corrections in g are due to 
the cut-off a), and we will get a solution similar to the known profiles of self-gravitating 
systems in the mean field approximation. On the other hand, if L ^> 1/mo, the spatial 
dependence of the solution will be washed out and we will have the reported constant 
solutions. However, the two solutions correspond to extremely high and extremely low 
densities, respectively. The system will only be homogeneous at such extreme densities. 
For intermediate densities, inhomogeneities will necessarily develop: the system will be in 
a mixture of high and low density domains. 

This instability as the volume increases had been notice long ago by Antonov, who 
found that the solutions of the Lane-Emden equation ceases to be (local) maxima of the 
entropy if the size of the box is larger than 0.335GAf 2 /(— E), where M is the total mass 
and E < the energy, and the system collapses |36j . This was called the gravothermal 
catastrophe in Ref. |3"7| . 
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6 Final remarks 



In thermal equilibrium, a self-gravitating system is made up of domains of low and high 
density (voids and clusters). The distribution of domains is a dynamical question that 
has to do with the way equilibrium is reached. The correlations within each domain are 
self-similar on a vast range of scales, decaying as 1/r. The transition to homogeneity 3 [3S] 
would take place on scales comparable to the range of the gravitational interaction, which 
may be larger than the deepest observed distances. A similar behavior (first order phase 
transition and self-similar correlations) was found for the two-dimensional self-gravitating 
system by using techniques of conformal field theory jSH] • 

It is straightforward to take into account the cosmic expansion in a simplified way, as 
in Refs. |14| I13j. by introducing comoving coordinates a: in Q}, such that the physical 
coordinates are r = R(t)x, where R(t) is the scale factor. As intuitively expected, this 
is equivalent to rescale the lattice spacing: a(t) = R(t)a and, consequently, we have the 
following rescaling of parameters: k — > R(t)n and ttTq — > R 2 (t)mQ. This variation of the 
parameters implies that the difference between the densities of the two phases decreases. 
This is not surprising: the expansion of the universe acts as a pressure that competes with 
the tendency of gravity towards collapse. 

It is remarkable that the picture of the self-gravitating system devised in this work 
is strikingly similar to the observed universe, despite it is not at thermal equilibrium. 
The scaling exponent of the the correlation function of the density, 7 = 1, however, is 
far from the widely accepted exponent of the galaxy correlation function, 7 ~ 1.8 for 
distances between 0.2 and 20 Mpc [Q. As mentioned in the introduction, in the last years 
there has been a strong debate about the validity of this result. Pietronero and his group 
analyzed the data from a different perspective and claimed that the exponent of the galaxy 
correlation function is 7 ~ 1, and that the self-similar behavior extends up to the deepest 
explored distances. The controversy seems not resolved, although most astrophysicist 
believe the earlier result, 7 ~ 1.8. This is theoretically supported by the thermodynamic 
arguments given in |T3]. These arguments, however, are somehow heuristic and their 
validity may be questioned. Indeed, using Eq. (|12|) as starting point, the same arguments 
can be applied step by step to the formulation of the self-gravitating system given in 
this work, leading to the same predictions of Saslaw's book. However, the results of this 
work, based on a rather rigorous treatment of the self-gravitating system, are completely 
different. 

Hence, the conclusion of this paper is clear: either Pietronero and coworkers are right, 
and 7 ~ 1, or the thermodynamical approach is not valid to describe the universe at such 
scales. The later possibility cannot be discarded, since it is difficult to argue that the 
universe is at thermal equilibrium, and dynamical effects may play a prominent role in 
the behavior of the galaxy correlation function. In any case, the results obtained in this 
paper need not be applied to the universe as a whole, or to the large scale structure of it, 
but to any self-gravitating system that are close to thermal equilibrium. The interstellar 
gas, where self-similar behavior has also been observed, might be an instance. In this case 
the scaling exponent of the density correlations is compatible with 7 = 1, but with large 
uncertainties [T5] . 

3 The scale at which the density fluctuations cease to be self-similar. 
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Appendix 

Let us outline in this appendix the derivation of equation 1)12(1 . which is not completely 
straightforward. The field <p n can be obtained from the derivative of the action as 

1 \ — i 1 \ -< — i dS 1/2 — h , 

1 E V nn>Pn> + ^ £ V^,— - -I , (20) 



K 

n' 



where we used Ylm 1 ^nri = V"^o- Now insert the above expression for <p n in equation (jTTj) 
The term proportional to (1/2 — K)j KfrL^ cancels the Ing term of Eq. (|llj) . Averaging over 
the thermal fluctuations the resulting expression for e, we are led to compute (p n dS / d(j) n i) , 
which can be written as 

Pn-Hr-) = -4 / [d^]Pn^-^ S ■ (21) 



4 



z „ 

Integrating by parts and taking into account that dp n /d(j) n i = 5 nn /p n (l — p n ) , we get 

dS \ 

Pn-^T~ ) = <W(/°n(l - Pn)} ■ (22) 
o<Pn' I 

The term linear in p n cancels the term proportional to k ^T) nn of the averaged Eq. (|1 ID . 
Hence, we obtain for the averaged energy 



e = k B T 



3 11 \ i \ If \ _i . \ 

^P ~ n 77 2^ K \Pn V nn>Pri) + 7j77 2^ K (Pn^nnPn) 



(23) 



nn' 

It suffices to write n' = n + r to realize that the above equation is Eq. (|12|) . 
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